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Abstract. 

The present work is devoted to the investigation of the interaction between 
vortices (topological defects) and site-impurities (structural defects) in the 2D XY 
model and its influence on the well-known properties of the pure system. The 
main goal is a theoretical description of the Berezinskii-Kosterlitz-Thouless (BKT) 
temperature reduction by quenched non-magnetic impurities, based on the vacancy- 
vortex interactions and the vortex-pair dissociation mechanism of the transition. The 
non-magnetic impurity interaction with a system of vortices can be found either from 
the phenomenological theory of topological defects or from the Villain model. We take 
both paths and compare the results obtained. Our prediction for the BKT temperature 
reduction is confirmed by the available Monte Carlo data. 
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1. Introduction 

An object of our interest is the two-dimensional XY model, described by the 
Hamiltonian: 

H = -jJ2 ( S r S ?> + S r S r>) ■ W 
<r,r'> 

Here, S r are unit spins placed on the sites of the square lattice with the nearest 
neighbour interaction J. This model is well known for its extraordinary properties 
connected with the presence of topological defects (vortices). It can also be 
considered as the limiting case of the 2D easy-plain Heisenberg model: i? e .p. = 
—JJ2<rr'> i^r^r' + S^S^i + XS^S^,), when A = 0. It has been argued that quasi-two- 
dimensional types of real magnetic materials, like layered magnets or ultrathin magnetic 
films can be satisfactory described within the 2D easy-plane Heisenberg model pQ, and 
since the behaviour of this model has been found qualitatively similar to that of the 2D 
XY model in a rather wide range of anisotropy parameter [21 E], it seems natural to use 
the 2D XY model as a suitable device for the study of real quasi-2D magnets. Thus 
the question of the influence of impurities, always present in the lattice structure of real 
materials, should be and has been posed in recent years [H EJ E] . 

We define here the 2D XY model without the non-interacting (and thus 
unimportant) component S* in the trace of the system, as is often done in literature. 
This case is also referred to by some authors as the planar rotator model (see, for 
example, [7j). This should not confuse the reader, since the non-interacting component 
would not change the qualitative picture anyway. 

The Hamiltonian ([1]) can be written in a more convenient form for calculation in 
terms of the angle variables —tt < 6 r < n: 

H = -JJ2J2 cos(9 r+aa -9 r ) (2) 

r a=x,y 

where & x = (a, 0) and = (0, a) form an elementary basis of the lattice, and J is 
the ferromagnetic coupling. Then, since we want to study a system with non-magnetic 
impurities (vacancies) in the lattice we introduce the "occupation numbers" : 

1, if there is a spin; . . 

0, if the site is empty. 

and construct the Hamiltonian: 

H = -J^2 Yl cos (^+a Q - # r )c r cv • (4) 

r a=x,y 

The introduction of such impurities (in a general case) makes the model impossible 
to diagonalize in the spin-wave approximation by a Fourier transformation as it is 
possible to do for the regular (without structural defects) lattice [8]. One distinguishes 
quenched and annealed types of dilution. Annealed dilution is understood as impurities 



Interplay of topological and structural defects in the 2D XY model 



3 



being in thermodynamical equilibrium with the spin degrees of freedom, so the averaging 
over the occupation numbers ([3]) should be taken already in the partition function [9]. 
In the quenched dilution case, the impurities are frozen at their positions with some 
fixed probability and one should average the observable quantities (like the spin-spin 
correlation function or the free energy of the system) over different configurations, 
and not the partition function itself. The latter statement, formulated in [9], was 
subsequently rigorously proven in [TO] . 

The random (quenched) dilution means that probability to remove a spin from a 
site is fixed and independent on the other sites state. So the averaging for all the possible 
configurations of vacancies can be written as 

~)= £ p (M)(-) > (5) 

{c r =0,l} 

with the probability function 

P({Cr}) = II[ C( W + (1-C)<W)] • (6) 
r 

This distribution is set in such a way that we obtain in average a system with 
concentration of magnetic sites c (or fraction of impurities 1 — c). 

Although the case of annealed impurities seems to be well studied and clear 
enough [HJ [12J, [13] , the influence of quenched dilution is a problem for which there are 
still unsolved questions and which deserves attention. For example, up to our knowledge, 
so far there are only Monte Carlo results for the phase diagram (Tbkt, c) j5][6], showing 
the reduction of Tbkt with decreasing concentration c of the magnetic component, 
and no theoretical constructions trying to explain this reduction. Also an approach to 
investigate the diluted model in the spin- wave approximation has been proposed in [14]. 

It is well known that the Berezinskii-Kosterlitz-Thouless transition is driven by 
the topological defects and is in some sense equivalent to the neutral 2d Coulomb gas 
transition to conducting state [15], [16]. In the present paper we describe the critical 
temperature reduction by an analysis based on the vacancy-vortex interactions (Section 
[3]). The form of this interaction can be found either from the phenomenological theory 
of topological defects [TT], [18] or directly from the Villain model [19] which can be 
regarded as a low temperature approximation of the 2d XY model [2TJ] (Sections [2] and 
H] respectively) . 

2. The vacancy interaction in a system with vortices 

An efficient way to study vortices in the 2D XY model is a continuous elastic medium 
approximation where the spin-wave excitations are forgotten and the topological defects 
are obtained from the "elastic" energy minimization under some special topological 
constraints. The spin variables 8 r defined on the sites of the initial lattice are promoted 
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to a continuous field 0(r), and the continuous limit of the spin-wave (harmonic) 
approximation of the Hamiltonian ([2]) is taken as the "elastic" energy of the system: 

E * = \ J H E ( A «^) 2 " \ J J rfr (W(r)) 2 . (7) 

r a=x,y 

The configuration of the field 9(r) that satisfies the topological condition § d9 = 2irq 
(definition of a vortex with winding number q), where the integral is over an arbitrary 
path enclosing the point defined as the vortex center, and has the minimal elastic energy 
(JTj), can be written in a polar coordinates system (centered at the center of the vortex) 
as: 

9 = qip + const , (8) 
q 

and its gradient, V9 = -(— simp, cos<^), can be found easily. This gradient is always 

r 

perpendicular to the radius- vector of the point drawn from the origin. The configuration 
obtained is called a vortex with the charge (winding number) q. The vortex is completely 
set by its charge, the constant in ([8]) is absolutely arbitrary, since one can switch from 
a configuration with one constant to a configuration with another constant without 
changing the energy (although the field configuration visually depends on the value of 
the constant). 

Actually, the total energy of such a configuration can not be correctly expressed by 
((7j), since in the continuous limit we have a singularity in the center of the vortex. Due 
to this, one has to specify the core energy of the vortex which is always finite and the 
elastic energy becomes: 

E^ Te = q 2 Jn / — = q 2 j7cln(L/A) . 

J A r 

It is divergent with the system size L and A is the radius of the core. We do not 
touch here the nontrivial question about the size of this core region and its energy 
estimation [18J. 

The elastic energy of a vortex with a non-magnetic vacancy at some sufficient 
distance r from the center can be found as the energy that corresponds to the four 
bonds removed (square lattice) subtracted from the energy of the pure system: 

£dil = ^pure _ E ^ = ^pure _ 1 j £ [(W . ^2 + ( _ W . ^2] (g) 

a=x,y 

= - \iJy 2 {2 sin 2 <p + 2 cos 2 if} = E^ c - Jq 2 (a/r) 2 . 

Thus, a non-magnetic vacancy has an attractive interaction with either a positive or a 
negative vortex charge. This is in good agreement with references [U [21]. Of course, 
this result is obtained via the assumption that the vacancy does not disturb the vortex 
configuration, an hypothesis which was reliably argued in [4j. Our result is almost 
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equivalent to that found in the paper mentioned, but seems to be a bit more definite 
since in [3] the coefficient in the interaction depends on the way of cutting out an area 
of the continuous field around the vacancy, and in our case it is only a matter of the 
lattice structure. 

We go further and consider a vortex pair (winding numbers q, q') containing an 
impurity. The spin configuration of the pair is simply the superposition of the single- 
vortex fields: 6(r) + d'(r) and thus the gradient is W + V0'. For the pure system the 
simple integration gives the elastic energy of such a pair: 

£P ure = -27T Jqq' ln(R/A) + ixJ(q + q'f \n{L/A) , (10) 

where R is the distance between the two vortices. Note that the second term, which is 
divergent, vanishes for a neutral pair (q' = —q). 

Let the polar coordinates of the impurity be (r, ip) in the coordinate system centered 
on the vortex q and (r', ip') in the system centered on the second vortex q' . We write 
down the result for the energy associated with this vacancy: 

£vac = Jo 2 ({q/rf + {q'/r'f + 2(q/r)(q' /r') cos(y? - if')) . (11) 

For a system with an arbitrary number of vortices and a vacancy at the point r we 
can generalize the elastic energy as: 



Ef = E*r - EW) (12) 

R R' 1 11 11 11 1 

where the sums span all the topological defects present in the system. 

So far we have been considering only one impurity in the system. Of course, a 
single vacancy does not have any influence on an infinite system, but, having formula 
(TTTj) for the vortex-pair-vacancy interaction, we can pass to the case of a finite fraction 
of empty sites and make some conclusions about the critical temperature behaviour with 
the concentration c. This will be the subject of the following section. 



3. The critical temperature reduction by quenched dilution 

In the previous section we obtained the form of the spinless site interaction with a 
system of topological defects. Instead of this, one considers now a single vortex pair in 
a system with some fraction of spins removed. The elastic energy can be written as 

^el^^r-E^acW, (13) 

*"vac 

where E^ re is the pure system energy, E vac (r) is the energy associated with a vacancy 
at the point r, given by ffTT]) . and the sum is over all the vacancies in the system. 

Obviously (flBT) is not exact because there are certainly some impurities which 
occupy neighbouring sites and thus destroy smaller numbers of links per vacancy. Their 
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contribution to the energy will be different, but (TT3"j) can be considered as a reasonable 
approximation when the concentration of dilution is weak enough. Eq. ffTB"]) can be 
equally written in the continuous approximation: 

Ef = E*r - J rfrp vac (r) J E; vac (r) , (14) 

with the impurity density introduced as 

Pvac(r) = ^5(r-r')(l-c r ) , (15) 

r' 

here 5 is for a delta-function and c r -s are the occupation numbers ([3]). 

The energy (TTJJ can serve to estimate the transition temperature, T BKT . Consider 
an ideal system that is constituted of a single neutral pair of vortices with winding 
numbers of modulus 1, thus one has only one degree of freedom, the separation 
R between the vortices One can define Tbkt as the temperature when this pair 
dissociates [22], i.e. the thermodynamical average 

, ,. f°° K'e-^VdR , , 

<* 2 > = <«9 

J a 

goes to infinity. With the undiluted system it happens at kT BKT /J ~ n/2, since 
E^ Te (R) = 27rJ\n(R/a) and one easily finds (R 2 ) = a 2 (nf3J - 1)/(tt/3J - 2). Using 
the best present estimate of the BKT temperature, kT BKT /J ~ 0.893 [23], it is 
obvious that this calculation gives quite a rough result. Nevertheless, in spite of 
a quantitative vagueness, this approach is qualitatively correctly based on the BKT 
transition mechanism. We choose to use it due to its simplicity and expect it 
to be satisfactory to examine the influence of non-magnetic dilution on the critical 
temperature. 

Now, with the energy (fl4"l) of a vortex-antivortex system with spin dilution, we can 
search the BKT point as the temperature where (fT6|) diverges. With p vac defined for 
an arbitrary configuration of impurities as (Tl5l) it is quite complicated, however one can 
use some approximate form of the density, reflecting its essential features. Here the 
quenched and annealed dilutions should be discriminated. In the frozen case the spins 
are removed randomly with the same probability for each site, so there is no preference 
for any part of the lattice to be more or less diluted than the rest of the system. Of 
course, fluctuations of random nature rather than thermal origin exist. The probability 
for these fluctuations goes to zero as the size of the lattice increases to infinity. Based on 
these arguments, a perturbation expansion has been proposed [H] where the 0th order 
can be considered as a "perfectly homogeneous" dilution. Here we get the corresponding 
approximation replacing (TT5|) with a "smeared" density: 

p(r) ~ (1 - c)N/(a 2 N) = (1 - c)/a 2 , 



which is simply the number of empty sites divided by the total volume. Now the integral 
in (}T4"|) can be simply calculated and gives E e \(R) = [1 — 2(1 — c)}27iJ\n(R/a). It follows 
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that the BKT temperature is just kT^ T /J = [1 — 2(1 — c)]tc/2, or, normalizing to the 
pure model critical temperature, 

^t/T b P kt = 1 - 2(1 - c) . (17) 

The critical temperature decreases with the dilution concentration, as one naturally 
expects due to the decrease of the average coordination number. Moreover, although our 
derivation was based on the assumption of a weak dilution, formula (17) predicts a 
vanishing of the Tbkt at c = 0.5. Being qualitatively correct, the last estimate differs 
from the known site percolation threshold concentration on a square lattice c ~ 0.59 




p = 1 - c 



Figure 1. The phase diagram of the 2D XY model with the quenched dilution of 
concentration p = 1 — c (c is the concentration of magnetic sites) . The Monte Carlo 
simulation results of [B] are compared with our theoretical prediction (1 1 7[) . The insert 
shows the vicinity of the percolation threshold. 

While the influence of quenched dilution in particular spin models, for example in 
the 2D and 3D Ising model [25, 26J, is well studied in numerous Monte Carlo simulations, 
for the model under consideration the computer experiment results are rather poor. 
We compare our result (1171) with the available Monte Carlo data for the diluted two- 
dimensional XY model [6] (Fig.l). The simulations were performed for the XY and 
planar rotator models with quenched dilution in two dimensions. Note that in terms 
of the paper [6j our model (HI) is the planar rotator model (PRM). Results of [6] for 
the PRM critical temperature estimated from two different methods, by the helicity 
modulus jump (PRM-Y) and by the spin correlation function exponent rj behaviour 
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(PRM-77), differ essentially in the region of weak dilution. The XY model points are 
again different. Although all three MC results go eventually to the correct percolation 
threshold (see the insert in Fig.l) there are few low concentration points and they do not 
seem to be reliable enough to make some strong conclusion about our result. However, 
at least the linear character of equation ([1711 seems to be present in all three MC sets 
in the weak dilution range, and this observation is also supported by [5]. 



4. The Villain model with nonmagnetic impurities 



According to [20] the Villain model can be derived in the low temperature limit from 
the Hamiltonian 



H 



-J I cos ( 

(r,r'> 



1 



;is) 



which is equivalent to that of the 2D XY model ([2]). We apply here the scheme of this 
derivation to the case of a diluted model, starting with the Hamiltonian 



<r,r'> 



cos(6> r — 9 r >] 



(19) 



with c r -s being the occupation numbers ([3]). The partition function of the model is then 
written: 



n 



n d9 r 



exp 



J2v(6 r - 9r>) c r c r 

<r,r'> 



(20) 



with V{6) = K [ cos 6 - 1 ] and K = J/k B T. 

The goal of the derivation is to find an approximate form of the expression under the 
integral in (T2"U|) . which preserves its initial symmetry and is easier to integrate. Namely it 
will be a superposition of exponents with quadratic arguments like in the SWA, but some 
new descrete variables will appear which subsequently can be associated with the vortex 
excitations in the system. Going through this procedure with the impurity variables c r -s 
one can expect to find the influence of dilution on the vortex energy contribution. 

As the first step one has to decompose the Boltzmann factor in (!20j) in Fourier 
series: 



cxp 



Y V(9 r - 9 r ,) c r c r 

(r,r'> 



n e e (*) e * 



S(dr-e r ,) „V(s)c r C r , 



(21) 



(r,r'> ■ 



where 6(s) = c r c r / + (1 — c r c r ')<5 Sj o naturally appears to ensure the equality when 
c r c r / = 0. The Fourier variable s depends on two real space variables: s = s(r, r'). 
Now, applying the Poisson summation formula [27] : 



00 



E 



m=— 00 



d(j) g{4>) e 



-2wi(j)m 



(22) 
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usually used to improve the convergence of the series, to each of the sums in (1211) . one 
can rewrite the partition function as: 



(lift) £ 



IZlr.r'S Vo(9r-0 T ,-2nm T y) c r c r , 

1 



(23) 



with e Vo ^ = J™ dcj) e v ^ e %< ^ 6 '. So far, no special assumption have been made and the 
result above is exact. 

Now let us consider the low temperature limit. In this approximation it is not 
difficult to find that e v °^ ~ e~ K02 ^. The latter comes as the result of the asymptotic 
form (K -> oo): e^ (s) = f ff" d6 e~ ise e K{ - cos() - l) « -JL- e ~ s2 ^ 2K \ One obtains the 
partition function 

Z= g e K. r ')(n/' r ^) e - K ^y)ier-e,-2, mrtVl ? CrCvl (24) 



of the desirable form, but as far as the limits of integration remain (— vr, tt) all terms 
with m^O give vanishing contribution (as K — ► oo) and (124j) is equivalent to the SWA. 
To repair this, extending the limits of integration, one can use the equality: 



d<pf(<p) = \im2^P^ / dip f(<p) , (25) 



true for any periodic function f(ip) — f(f + 2n). Eq. (l25l) can be easily checked by 
passing to the Fourier transform: f(cp) = Yl^L-oo e iS!p F(s). The left part of (12"5T) is just 
27rF(0). Integrating term by term the right part of fl25|) and taking the limit e — > one 
finds again 2irF(0). 

Finally, one has the partition function that describes the Villain model with non- 
magnetic impurities: 

Z= £ 6K, r ,) (J! r ^) e-^v/n, (26 ) 

m(r,r')=-oo V r " / -°° / 

with the Hamiltonian 

H™ } = J (fir ~ 6 r' - 2lTm r y) 2 C r C r/ + ej^^r ' ( 27 ) 

<r,r'} r 

When all the c r -s are taken equal 1 it turns to the Hamiltonian of the regular Villain 
model [IS]. 

It is known that the pure Villain model Hamiltonian can be divided into two parts: 
one corresponding to the spin-wave excitations and another one that corresponds to the 
vortex contribution. It is achieved by the Fourier transformation of the spin variables, 

* = ^E *. = ^E^. 

v k v r 

Fourier transformation of the discrete variables m rr >, 
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m r , r+aa = J= e- Jq(r+a ^ /2) m^ , a = x, y 
q 



— V e iq( 



I' 

and the shift of the Fourier transform of the spin variable: 



$k — Pit — 27TZ- 

with K a (k) = 2sin^f. 

Applying this scheme to the diluted Hamiltonian ( J27l) we find that again, as in 
the pure case, the two parts - the spin-wave and the vortex one - can be distinguished 
but now a third term appears which depends both on the spin and vortex degrees of 
freedom. Thus the Hamiltonian is made of three terms as: 

Zjdil , rrdil , rxdil {OQ\ 

-"Vill ~~ -"SW -"Vort "r -"SW^ort ■ 

Of course, the cross-term, Vort , vanishes in the pure model limit: c r = 1, r = 1, N. 
It has the form: 



#sw,Vo rt = 4vr J P( k + k (^( k + k')^(k)^(k') 

k,k' 

- L y (k + k')^(k)^(k')) (^ 2 (k') + ^ 2 (k')) _1 ^k' , (29) 

where = i{Ky(k)m^ — /^(k)??^) are the Fourier transforms of the vortex charges, 
L Q (k)=cos*f, K a (k) 
was defind before, and 

^) = ^E eJqr ( 1 - C r)- ( 3 °) 
r 

The parameter ( 1301) characterizes the strength of dilution |5]. When there are no 
impurities in the system one gets p(k + k) = and thus -ffgw,Vort = 0- The spin- 
wave and vortex parts of the Hamiltonian also contain terms which depend on p and 
which turn to zero in the pure case as well. 
The spin-wave term: 

k a 

+ J^p(k + k') (^L Q {k + k')K a {k)K a {k') ) cfrtpv , (31) 

k,k' \ a J 
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contains the "pure" spin-wave Hamiltonian (first term) and a contribution of the dilution 
that vanishes when p(k + k') = 0. The dilution contribution naturally has exactly the 
same form as was found in [5]. 

The new result here is the form of the vortex energy in the presence of non-magnetic 
dilution: 



< rt = 2vr J J k L)^ + 47r J E ^ k + k ') ( 32 ) 



k,k' 



L x (k + k')i^(k)i^(k') + L y {k + k')K x (k)K x (k') 



<?k<?k' 



(^ 2 (k) + ^ 2 (k))(^( k0 + ^(k')) 

Again one has the vortex interactions similar to those of the pure Villain model (first 
term) while the dilution effect is represented by the second term. It is known for the 
"pure" Villain term that: 



-2 T yj«(R)«(R')ln(|R-R'|/a) + ir 2 J^( 9 (R)) 2 , (33) 



R,R' R 



where R, R' span the sites of the dual lattice and g(R)-s are the vortex charges or 
winding numbers [T9] . 

Note, that in order to present the dilution contributions of the Hamiltonian (!28j) 
in an easily readable form we have not included into Eq. ([29l) . (13T1) and ( 1321) the terms 
quadratic in p. Anyway, one can neglect them in the approximation of a weak dilution 
(see [5]) which is the case here. 

The expression presented in ( 1291) and especially the form of the impurity 
contribution to the vortex part of the Hamiltonian, Eq. (!32l) . can further serve to analyze 
the impact of dilution on the peculiarities of the BKT transition. Let us first find an 
approximation that would correspond to the "smeared" impurity density approximation 
of Section 3. 

Imagine that the fraction 1 — c of sites is removed in such a way that the vacancies 
form some regular structure, then 

p(k + k') = (l-c)5 k+ k',o- (34) 



Of course, with random dilution this is not the case, but, as was argued in [T4|, (J3] 
can be considered as the zero approximation when one neglects the random fluctuations 
(homogeneous dilution). Moreover, (1341) is the disorder averaged value of p. We will 
see that this replacement with its average value corresponds to the "smeared" impurity 
density approximation of Section 3. In this case we obtain the vortex energy: 

H™ t = 2tt [1 - 2(1 -c)]jJ2 v^fS^ • ( 35 ) 

k ^ 2^ 7 A 7 W 



Interplay of topological and structural defects in the 2D XY model 



12 



As a consequence the energy of a neutral vortex pair is E^(R) = [1 — 2(1 — 
c)]2nJ\ia(R/a), exactly the same result as what was found from the topological defect 
theory approach under the assumption of the "smeared" density of vacancies. This leads 
of course to the same estimate of the critical temperature as well. 

5. Conclusions 

We have exploited two different approaches to account for the influence of quenched 
impurities on the vortices in the two-dimensional XY model: in the frame of the 
topological defects phenomenological theory and from the Villain model Hamiltonian. 
The interaction of a vacancy with vortices in the theory of topological defects was found 
to be attractive, in good accordance with other works on this subject jH [21]. In order 
to estimate the critical temperature change we used an approach based on the vortex- 
pair dissociation mechanism of the BKT transition. The "smeared" impurity density 
approximation leads to the same predictions for the critical temperature within the two 
approaches (topological defects theory and the Villain model). The dependence of the 
transition temperature on the magnetic sites concentration, Tgj[ T (c), obtained under 
the assumption of a weak dilution, however gives a percolation threshold which differs 
from the known real site percolation threshold for a 2D square lattice. Comparing with 
the currently available Monte Carlo results |6J which unfortunately are not accurate 
enough to make reliable conclusion about the weak dilution range, we recover the linear 
character of the critical temperature decrease close to c = 1. 
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